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The purpose of this paper is to present general relativistic cosmological hydrodynamic equations 
in Newtonian-like forms using the post-Newtonian (PN) method. The PN approximation, based on 
the assumptions of weak gravitational fields and slow motions, provides a way to estimate general 
relativistic effects in the fully nonlinear evolution stage of the large-scale cosmic structures. We 
extend Chandrasekhar's first order PN (1PN) hydrodynamics based on the Minkowski background to 
the Robertson- Walker background. We assume the presence of Friedmann's cosmological spacetime 
as a background. In the background we include the three-space curvature, the cosmological constant 
and general pressure; we show that our 1PN approach is successful only for spatially flat cosmological 
background. In the Newtonian order and 1PN order we include general pressure, stress, and flux. 
The Newtonian hydrodynamic equations appear naturally in the OPN order. The spatial gauge 
degree of freedom is fixed in a unique manner and the basic equations are arranged without taking 
the temporal gauge condition. In this way we can conveniently try alternative temporal gauge 
conditions. We investigate a number of temporal gauge conditions under which all the remaining 
variables are equivalently gauge-invariant. Our aim is to present the fully nonlinear 1PN equations 
in a form suitable for implementation in conventional Newtonian hydrodynamic simulations with 
minimal extensions. The 1PN terms can be considered as relativistic corrections added to the well 
known Newtonian equations. The proper arrangement of the variables and equations in combination 
with suitable gauge conditions would allow the possible future 1PN cosmological simulations to 
become more tractable. Our equations and gauges are arranged for that purpose. We suggest ways 
of controlling the numerical accuracy. The typical 1PN order terms are about 10" 6 ~ 10" 4 times 
smaller than the Newtonian terms. However, we cannot rule out possible presence of cumulative 
effects due to the time-delayed propagation of the relativistic gravitational field with finite speed, 
in contrast to the Newtonian case where changes in the gravitational field are felt instantaneously. 
The quantitative estimation of such effects is left for future numerical simulations. 
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I. INTRODUCTION 



The nonlinear evolution of large-scale cosmic structures is usually investigated within the framework of Newtonian 
theory either in analytic studies or numerical simulations. Without investigating general relativistic effects, however, 
it is not clear whether the Newtonian theory is sufficient to handle such large-scale cosmological structures. If we 
regard Einstein's gravity theory to be the correct framework on such cosmic scales, the relativistic effects should 
exist always. The point is to which level we can practically ignore relativistic correction terms considering currently 
available levels of observations and experiments. In case relativistic or nonlinear effects are not large, we have two 
widely known ways to estimate relativistic effects: one is the perturbation approach [1-3], and the other one is the 
post-Newtonian (PN) approximation [4-9]. 

In the perturbation study, we recently have investigated the weakly nonlinear regimes based on Einstein's gravity. 
We showed that, except for the gravitational wave coupling, the relativistic scalar- type perturbations coincide exactly 
with the Newtonian ones up to the second order [10]. The pure relativistic correction terms appear in the third order. 
These terms turn out to be independent of the horizon scale and small (~ 10~ 5 order) compared with the second-order 
Newtoninan/rclativistic terms [11]. Our studies have shown, from the view point of Einstein's gravity, that Newtonian 
gravity is practically reliable near the horizon scale where structures are supposed to be weakly nonlinear. Therefore, 
justifying its use in current large-scale numerical simulations which cover the Hubble volume. Considering the action- 
at-a-distance nature of Newton's gravity our result is a rather surprising one because the relativistic perturbation 
theory is applicable on all scales including the super-horizon scale. Thus, our perturbation study ensures that to the 
weakly nonlinear order Newtonian theory is reliable even near and beyond the horizon scale. 

Now, how about situation in the nonlinear regime on scales smaller than the current horizon? On such scales the 
structures could be in a fully nonlinear stage. However, if the relativistic gravity effect is small we can apply another 
approximation scheme which is well developed to handle isolated bodies in Einstein's theory: the PN approximation. 
The PN approximation has been important to test Einstein's gravity. It presents the relativistic equations in a form 
similar to the well known Newtonian equations. In the PN method, by assuming that the relativistic effects are 
small, we expand relativistic corrections in powers of v/c. In nearly virialized systems we have GM / (Rc 2 ) ~ (v/c) 2 . 
Thus, the PN approximation is applicable in the slow motion and weak gravitational field regime. The first order 
PN (1PN) approximation corresponds to adding relativistic correction terms of the order (v/c) 2 to the Newtonian 
order terms. In this approach we recover the well known Newtonian hydrodynamics equations in the OPN order 
[5]. The PN approximation is suitable to study systems in which Newtonian gravity has a dominant role, while the 
relativistic effects are small but non-negligible. Well known applications include the precession of Mercury's perihelion, 
and various solar system tests of Einstein's gravity theory like the light deflection [12]. Recent successful applications 
include the generation of gravitational waves from compact binary objects, and the weakly relativistic evolution stages 
of isolated systems of celestial bodies [13]. Another important application is relativistic celestial mechanics which is 
required by recent technology driven precise measurements of the solar system bodies [14,15]. Notice that all these 
applications are based on the PN approximation of isolated systems assuming the Minkowski background spacctimc. 

In this work on the cosmological PN approximation we assume the presence of a Robertson- Walker background. 
Thus, we have the Friedmann world model built in as the background spacetime. We will show that the Newtonian 
hydrodynamic equations come out naturally in the OPN order which is the Newtonian limit. In the context of large- 
scale cosmic structures PN effects could be essentially important in handling gravitational waves and gravitational 
lensing effects. It is well known that Newtonian order treatment of the gravitational lcnsing shows only half of the 
result from Einstein's theory which is due to ignoring the 1PN correction in the metric. We also anticipate that 
PN correction terms could affect the dynamic evolution of large-scale structures, especially considering the action- 
at-a-distance nature of the Newtonian theory. This can be compared with the relativistic situation in which the 
gravitational field propagates with finite speed. In this work we present the complete set of equations to study 1PN 
effects in a cosmological context. The PN approach can be compared with our previous studies of the weakly nonlinear 
regime based on the perturbation approach which is applicable in the fully relativistic regime and on all scales. In 
comparison, although the PN equations are applicable in weak gravity regions inside the horizon, these are applicable 
in fully nonlinear situation. Thus, the two approaches are complementary in enhancing our understanding of the 
relativistic evolutionary aspects of the large-scale structure in the universe. 

In the PN approximation we attempt to study the equations of motion and the field equations in Einstein's gravity 
in a Newtonian way as closely as possible. The Newtonian and post-Newtonian equations of motion follow from 
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the energy-momentum conservation. The Newtonian order potential and 1PN order metric variables are determined 
by Einstein's equations in terms of the Newtonian matter and potential variables. Thus, the metric (or relativistic) 
contributions are reinterpreted as small correction terms to the well known Newtonian hydrodynamic equations. The 
Newtonian equations naturally follow in the OPN order. The form of 1PN equations is affected by our choice of the 
gauge conditions. In this work we take unique spatial gauge conditions which fix the spatial gauge modes completely; 
under these spatial gauge conditions the remaining variables are all equivalcntly spatially gauge- invariant, see Sec. 
V. The temporal gauge condition (slicing condition) will be deployed to handle the mathematical treatment of the 
equations conveniently. We show how to choose the temporal gauge condition which also removes the temporal gauge 
mode completely. Under each of such temporal gauge conditions the remaining variables are equivalently (spatially 
and temporally) gauge-invariant. We arrange the equations and gauges so that the final form of equations is suitable 
for numerical implementation in conventional cosmological hydrodynamic simulations. 

We present basic quantities and the derivation of our 1PN equations in Sec. II and III. If the reader is more 
interested in the application of the 1PN equations, the basic set of cosmological 1PN equations is summarized in Sec. 
IV. The gauge issue is expounded in Sec. V. We discuss our results and several unresolved issues in Sec. VI. 



II. BASIC QUANTITIES 
A. Curvature 

As the metric we take 



5oo = - 



1 - \ (2U 2 -4$) 



o- 



9H = a 2 (l + ^TVj-Yij +0~\ (1) 

where x° = ct, and a(t) is the cosmic scale factor of the background Friedmann world model. Indices a, b, c, . . . 
indicate spacetime, and k, . . . indicate space. Tildes indicate spacetime covariant quantities, i.e., spacetime indices 
of quantities with tilde are raised and lowered with the spacetime metric g a b- The spatial index of Pi is based on 
jij in raising and lowering indices with 7 y , an inverse of jij. jij is the comoving (time-independent) spatial part 
of the Robertson- Walker metric, see Eq. (2) in [16]. In a flat Robertson- Walker background jij becomes % if we 
take Cartesian coordinates. Compared with our previous notations used in perturbation studies in [17], we set the 
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comoving part of the three-space background metric g\- = jij and use part of the Latin indices to indicate the space. 
We are following Chandrasekhar and Nutku's notation [5,7] extended to the cosmological situation, see also Fock [4]. 
The 2U/c 2 term in goo gives the Newtonian limit, and if we ignore all the Newtonian and post-Newtonian correction 
terms we have the Robertson- Walker spacetime. Thus, our PN formulation is built on the cosmological background 
spacetime. 0~ n indicates (v/c)~ n and higher order terms that we ignore. The expansion in Eq. (1) is valid to 1PN 
order [5]. Dimensions are as follows 

[g ab } = [9 ab ] = l, [li j ] = h ij ] = ±, [a] = l, [c]=LT~\ [U] = [V] = [c 2 ], [PJ = [f*] = [c 3 ], m = [c% (2) 

where L and T indicate the length and the time dimensions, respectively. 

In our metric convention in Eq. (1) we have ignored the possible presence of ^ + Cjlj + Cjli) like terms in 

(jij with C 1 ^ = by choosing the spatial C-gauge conditions (C = = Cj) which remove the spatial gauge mode 
completely; this will be explained in Sec. V; a vertical bar indicates the covariant derivative based on jij. Under such 
spatial gauge conditions, we can regard all our remaining PN variables as equivalently spatially gauge-invariant ones, 
see Sec. V. We still have a freedom to take the temporal gauge condition which will be chosen later depending on 
the mathematical simplification or feasibility of physical interpretation of the problem under consideration, see Sec. 
V: in the perturbation approach we termed this a gauge-ready approach [18]. We also have ignored the spatially 
tracefree-transverse part of the metric (CV, with C| = = C^ ) because gravitational waves are known to show up 
only from the 2.5PN order [8]. 

The inverse metric becomes 
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-00 = 



~0i 



The determinant of the metric tensor g is 

where 7 is the determinant of 7^ . 
The connection is 



1 + i 2C/ +i( 2 ^ 2 + 4<&) 



3 n-^^f+o- 4 



f° 

1 00 



l + ^(3F-i7) + 0- 



(3) 



(4) 



-2$ + -P^C/i + iT 1 <ir 7 , 
a 



fo 

1 0i 



\u ti - 1 (2$,, + aPi) + L-'O- 6 , 



a 1 



1 00 — 



tU' 



1 a 



1 Qj — o,- i 5- 



v + 2^ ([/ + v) J 7 « + -p m \ I + l- 1 ©- 5 , 

2(U + V) V* - 2$'* - (aP 1 )' J + L~ 1 (D~ 6 , 



+ z,- 1 er 5 , 



rj* = r ( ^; fe + 1 (v fe ^ + « - v^ jk ) + l-'o- 



(5) 



r' 7 ^- /i , is the connection based on 7^; we introduce P(i|j) = \ (Pi\j + Pj\i) and P[%\j\ = | (Pj|j — Pj|i)- We have 
1 du — if T 

c dt — c u ■ 

The Riemann curvature is 



where T 

JJ _ 1 du - 1 JT 



R° 0Qi = -1 (aPi + -U Aj P^j + L- 2 0-\ 



R° Qii = L-'O- 6 , 



R° — 
11 iOj — 



+ u Aj ) 



V + 2-(U + V) 



la + 



-aaU + 2d 2 (U + V) 



-2U t{i V tj) - UiUj + U' k V k7l3 + 2* <y + {aP m y j + L^O' 6 , 



R a 



ijk 



R l OOj ~ 



2[V+-U 

a 



1 2K 

lk]i - - P [j\k]i + —li\j P k] 



+ L~ 2 0~ 5 , 



1 / a „■ 1 rr j 



1 



V/ + ^(f7 + 2y) 



<f- 

J a 2 



+ 



2a 2 



2(C/ + T/)C/' l b . 



+ L~ 2 



rt Ojfc - -3 



2 5t,--P, 



+ z,- 2 e>- 5 , 



Lfi _ 



M - ^S-* + 1 M u^ jk ) + 1 (p*„ - p, '<) |fc 

& jkl = -2^ 7j[fc <5i ] + 12 (-a 2 ^ + V m 8\ V>\ [kll]3 ) + L- 2 0-\ 



+ L- a O 



2/n-5 



(6) 
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It is convenient to have 



pi — pi r(t)* p' P — P _l p(T)' p 

^ \jk = ^ \kj U Ijk^ ' Mb'fc - ^i\kj + n t jk r l> 



R 



jkl 



K{8i ljl -6\ ljk ), R^=2K 7ij , R™=6K, 



(7) 



where K indicates the comoving (time-independent) part of background spatial curvature with dimension [K] = L 2 . 
The Ricci curvature and the scalar curvature become 

R° = 1 ( 3 | + ±U) + 1 {3V + 3^ (f/ + 2V) +6^-1 (C/ - V), + 2VAU - 2A<I> - (aP^ 

+L- 2 0- 6 , 



2 ( V+-U 

a 



— ,P j :..-AP i -2KP i ) 
2a V V 



R'n ~ T ;t 



2 y 



>a ( 



L- 2 0-\ 



P\ 3 l - AP i + 2KP i 



R) 



R = 



a 2 3 c 2 
6K 1 



a a 

- + 2 - 

a a, A 

• 2 

a a 

- + — 

a a- 2 



A + 4K\ 



a" 
A . 



+ L- 2 Q-\ 



+2—U-A 



A + 3K 



V 



+ L- 2 0~\ 



where A is a Laplacian operator based on jij . 
The Weyl curvature is introduced as 



C\ cd = R\ cd - \ (5 a c R bd - 6 a d R bc + g bd R a c - ~g bc R a d ) + ±R (6 a c g bd - 5 a d ~g bc 



To 1PN order we have 
1 K 



° OOi — 



C° 



Oij 



---P l+ L- 2 0-\ 
c J a 

L- 2 Q-\ 



1 



iOj 



\{U + V) Aj -\A{U + V) lij 



+ L~ 2 0- 



C" 



C 



1 a 
c"3 2 
1 1 



(p< 



IP 



AP [fc - 2XP r , 



[fe I 



2B 



OOj 



Ojk 



\(U + vf ] - l -A{U + V)5) 



° jOk - 



^ jkl 



1 1 

c 3 " ~2a 
1 1 

c 3 4a 



P' 



|i[fe -AP [fe + 2^P [fe )4-2P b . |fe] 
Af* + 2*^)7^- (P\ h - 



[j\k]i +L 2 °, 
+ L- 2 0-\ 

+ L- 2 Q- 5 , 

AP j + 2KP 3 ) 5\ + 2 ( 





/ Ife 



+ L~ 2 0- 



1 



-A ([/ + V) 5\ kll]j + (U + V) Mk 5\ - (U + vy i[k 



+ L- 2 0-\ 



(8) 



(9) 



(10) 



We can check the tracefree nature of the Weyl tensor: C b bcd = = C c bcd - Here we encounter non-vanishing traces 
involving K terms, like 



bOi 



±*p l+ L- 2 0-\ & 
c J a 



Oci 



±™ Pt + L- 2 0-\ 
c J a 



whereas 



& Qc0 = L- 2 O-\ C b bij =L- 2 0~ 4 , C c lc0 = L- 2 O-\ C\ c ^L- 2 Q-\ 



(11) 



(12) 



This indicates that the presence of K terms in the 1PN order appears to be not reliable. Later this point will be 
resolved as we show that the K terms indeed can be related to the higher order terms in the PN expansion, see below 
Eq. (88). Until we reach such a conclusion about the effect of K term we will keep it in our equations. 
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B. Energy-momentum tensor 



The normalized fluid four- vector u a with u a u a = — 1 gives, to 1PN order, 



u 



1 + 



1 /l 



; 1 ;+^r 



V + « 2 



F + F ] +^F 2 + 2$-F^ 



H i-o 
c a 



+ ^ + V +-F 



1. 



2$ 



+ er 6 , 







1 





















Pi 



o- 



(13) 



where the index of v l is based on 7^-; corresponds to the peculiar velocity field. The energy-momentum tensor is 
decomposed into fluid quantities as follows 



Tab = QC Z ( 1 + ft ) u a u h + p (u a u b 



-gab) + q a u b + %u a + n ab . 



(14) 



The energy density fi{= gc 2 + gfl) is decomposed into the material energy density gc 2 and the internal energy density 
gfl; p, q a , and TT ab are the isotropic pressure, the flux, and the anisotropic stress, respectively. We have q a u a = 0, 
TT ab u b = 0, 7T^ = 0, and TT ab = n ba , thus 



We introduce 



1 1 „ i 1 1 . , , 1 1. h 

90 = QiV , 7T i = TTijV J , 7T = o K^ij ' 

ca ca c z a z 



qi — -aQi, jtij = a 2 Hij, 



(15) 



(16) 



where indices of Qi and ILj are based on 7^ . We take q a and ft a b to have post-Newtonian orders as introduced in Eq. 
(16), see for example [19]. This will be justified by the energy conservation equation in the Newtonian context which 
will be derived later, see Eq. (102). In a strictly single component situation we can always remove Qi by following the 
fluid element. However, there exist situations where we have the additional flux terms present even when we follow 
the fluid elements; the fundamental origin of such flux terms can be traced to the presence of additional components 
in the energy-momentum. Thus, in our case it is more general and convenient to keep the flux terms separately. Up 
to the 1PN order the condition 7r£ = gives 



n* = ^ifff. 



(17) 



We also set 



g = g, II = IT, p = p. 



(18) 



Dimensions are as follows 



[«„] = [u a ] = 1, [fab] = [f b a ] = [f ab ] = [p] = [q a ] = [q a ] = [n ab ] = {gc 2 } = ML-'r- 2 , 

N = K] = [c], [n] = [c 2 ], [Qi] = m = [ cqa ] = [ Q c% [it,] = [n;.] = TF] = [n ab ] = [gc 2 ]. 



We have 



(19) 



Too = Qc 2 { 1 + - (v 2 - 2F + n) + - 



-- (2Q i v i + Ilij-vV) 



c 4 

Q-6 



v 2 [v 2 + 2V + - 2Fn + 2F 2 - 4$ 
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,1 1 

T 0i = -age < -Vi + — 



v 2 + 2v + n + v - - Pi + - (Qi + n ijV j ) 



+ o- 



fij = a 2 gc 2 \ 1 (viVj + ^7ij + ^Ily 



+ 



viVj (v 2 + 2u + w + n + Q + 2V^ 7ij - 2v {l p j) + -^Qavj) 



T = -gc z 

From Eq. (14) covariantly we have f = -gc 2 ^1 + -^flj + 3p. 

C. Fluid- frame kinematic quantities 

The projection tensor h a b based on the fluid-frame four-vector u a is defined as 

Kb = g a b + u a u b . 

The kinematic quantities are defined as [20] 

1 Z7 



(20) 



ab — h c a h1u c - d , 9 — U a . a , a a b = 9( ab } - -8h ab , U a b = 6[ab]i a a — U a = U a - b U b , 



(21) 



(22) 



where we have u a 9 a b = 0, u a a a b = 0, u a Lo a b = 0, u a a a = 0, and 9 = 9 a a . The kinematic quantities 9, a a , d a b, and ui a b 
are the expansion scalar, the acceleration vector, the shear tensor, and the vorticity tensor, respectively. 
To 1PN order, the projection tensor becomes 

h° = -^v 2 - i [v 2 (v 2 + 2U + 2V) - v'f*] + Q-\ 

h\ = a S^Vi + 1 [ Vi (v 2 + 2U + 2V) - Pj J + 0~ 5 , 

[v 2 (v 2 + 2U + 2v) - tfPj] | + er 7 , 



/in = U <^ 1 + ^ 



c a 



c 4 



fcj = 5) + + 1 (« 2 + 2C/ + 2V) - + O- 6 . 



(23) 



The kinematic quantities become 



ia (v^ + ajij) + ( a «i)' + Q" 2 + u + 2vj + (viV k \j + VjV i][k ) 

Jij +2v [i (U + V) tj] -P m \+L- 1 0- 5 , 



V +±{U + 2V) + -v k V k + l -*v 2 
a a 2 a 







1 















w + 3 ^ + 2 l ,,„. + I w ., + _L (a v )+ i- („*„<)„ 



I/O- 5 



(24) 
(25) 



ViVj [a+^-v' 



(v 2 v k ) 



+ L- 1 0- r \ 



-a(v 2 y + {U + 2V)v\ k + - 
u^ ] = \av m + la |«y (a« t] )' + Qv 2 + C/ + 2V^j + v k (v kl[jVi] + v [jVi]]k ) + 2v {l (U + V) ^ - P m j 



(26) 



(27) 
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a.i = — \lJ.i I 1 



1 



1 



-2$ 



cr ' \ c 
+L- 1 C- 6 . 

From u c a c = and u b 9 a b = we have 

1 1 



1 fl 



a 



— v l a i7 8 0i 
c a 



-v z + U 



1 1 



c a 



_a i 

at 



i i 



J%3 , 



7 00 — ~n 



2 1 VV ^ih 



and similarly for d$ a and ojo a . 

The electric and the magnetic parts of the Weyl (conformal) tensor are introduced as 



E ab = C acbd u c u d , H ab = ^h ac ef C, 



efbdU C U d , 



(28) 
(29) 

(30) 



1 ^abcd 



i "* is the totally antisymmetric tensor density with e d , the totally antisymmetric Levi-Civita 



where jj abcd 

1 V-g 

symbol with e 0123 = 1. E a b and H a b are symmetric, tracefree and u a E a b = = u a H a b. To 1PN order we have 



El 



-C 



00j 

1 1 



1 M 1 



(u + vyr--A(u + v)5* 



u°u° 



1 

~2' 
1 1 



f r c 



1 1 



kOj 



v l C° 



c a 



kOj 



+ L- 2 0-\ 
V lk ca VC 



OOj 



ikl 



\ (P m vnl AP l + 2KP t ) + ± Wl A (U + V) lk3 + P m -v t (U+ V) w } + L- 2 



(31) 



(32) 



where we introduced rfi k = with e yfe , the totally antisymmetric Levi-Civita symbol with e 123 = 1. Thus 

fjOijk _ ^J^yjijk _ ixidices of rfi k are based on 7^. E 0a follows from u b E a b = which gives E 0a = — (t? /u°)E ia = 

— (v l /ca)E ia , thus E ai ~ L~ 2 0~ 3 and E m ~ L~ 2 0~ 4 . For nonvanishing if, £r a6 is not tracefree to L~ 2 0~ 3 , but as 
we mentioned earlier the term involving K is already 0~ 2 order higher, thus consistent. 



D. Normal-frame kinematic quantities 

The normalized normal- frame four- vector n a , which is normal to the space-like hypersurfaces, is defined as hi = 
with h a h a = — 1. To 1PN order we have 



h° = l + ^rU 



1 



1 fl 



-U 2 + 2<S>)+0 



U 2 - 2$ + O 



1 1 



hi = 0. 



c 4 V2' 

By setting hi = in Eq. (13) we also recover the normal vector, thus 

Vi = \p i + LT- 1 0-\ 

where Vi is, say, the velocity of the normal frame vector. The projection tensor based on h a becomes 

hii=9n, h 0i =goi, h o = 0- 6 , ft*=5j, fc« = - + , ft? = = 4 
Kinematic quantities based on n a become 

6 = -3- + 1 f 31/ + 3-17 + -pO + i- 1 ©- 5 , 
c a c J \ a a ' J 



LUij = 0, 

hi 



c z c 4 



(33) 

(34) 
(35) 

(36) 

(37) 
(38) 
(39) 
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(Toe follows from a ac h c = 0, thus a M ~ L _1 C -6 and <t o ~ L -1 © -9 ; we have u a b — 0. Similarly, do follows from 
a c h c = 0, thus 5 — L _1 -5 . 

The electric and the magnetic parts of Weyl tensor based on h a give 



-& 



OOj 



1 1 



+ L~ 2 0~ 



- } n n r] C jH 



1 1 



/7.7 



APi + 2KP l ) l3k + P m 



\ml " 1 " i 



-2/n-5 



L- z O 



(40) 



(41) 



Thus, i?j is the same in both frames to the 1PN order; compared with Eq. (32) i?j differs. The nonvanishing K term 
which causes H a b to be not tracefree to L~ 2 (D~ 3 can be ignored because the K term is of the O" 2 order. 

E. ADM quantities 

In the Arnowitt-Deser-Misner (ADM) notation the metric and fluid quantities are [21] 



.goo = -N 2 + N l N t , g oi = N,, ~ 9lj = h ih 
ho = -N, hi = 0, h° = N~ 1 , h'^-N^N 1 

E = n a n h f a \ Ji = -n b Tf, S tJ =f tJ , S = h ZJ 5, 



Sij — Sij ^hijSj 



(42) 
(43) 

(44) 



where TVj, Ji and Sij are based on hij as the metric, and is the inverse metric of hij. The extrinsic curvature is 

1 



hij,o) , K = h %3 Kij, Kij = Kij — -hijK, 



(45) 



where indices of Ka are based on h^; a colon ':' denotes the covariant derivative based on ha with T 



(h)i 



\h tl {hj Lk + ha-j — hj k _i). The intrinsic curvatures are based on the metric h^ 

r>(h)i _ -p(h)i _ -p(h)i -p(h)m-p(h)i _ p(/i)m -p(h)i 

n jkl — 1 jl,k 1 jk,l +L jl L km 1 jk L Im' 



R (h) = R (h)k 
n i 3 - n ikj ' 



RW = WR$\ R^^R^-lhijR^. 



A complete set of ADM equations will be presented in Sec. Ill C. 
To the 1PN order we have 



jk 



(46) 



AT = 1 - 1[/ + 1 (hj 2 - 2$j + 0-\ 

N i = ~±P i (l-\u)+0-\ Ni = -\aPi 



1 + \(2V-U) 



h tJ = a 2 1 + -2V + 0-\ h^ 



-,(/>)» 

jk 



1 



2V + O- 



R (h) = R h) 



(V i] j+jijAV)+L- 2 0- 



R (h) 



6K - — 4 (A + 3K) V 



Ki 



K 



V+-(U + 2V) 



-y ij + -P m \+L- 1 0- B , 



l_ 3 a l_ 

c a & 



;-;( • . -f) . -I- 



+ L- 1 0-\ 



L- 2 0-\ 



(47) 
(48) 
(49) 

(50) 



9 



E = gc 2 i + 4^ + n) + l 



v 2 [v 2 + 2U + 2V + n + " ) + 2UU - 2Py + - (2Q i v i + U ij v i v j 



Ji = age 



f 1 1 

[c C 



r \ v 2 + u + 2V + n + Q - Pi + i (Q, + n y v) 



+ 2U + 4 V + n + i- + 2i-F 7ij - 2v {i P j) + -Q^Vj) 



S=QC 2 



a 2 (v 2 + 2U + 2V + n + Q - 2iV + i (2Q i w < + n«t;V) + J 



(51) 



III. DERIVATIONS 



A. Equations of motion 



Using Eq. (20) the energy and momentum conservation equations give 
1^» 



o = --fi b 



1 , 3 v 1 

— a 3 <x +- 
a 01 a 



1 



(TV 



1 

+ -^g 



V + -v l (V -U).+ V - t 

a '* a g 



+ e T- l o- 



0= -T* 

a *' 

1 



CTOj + — (Qi + n^V) 



+ - | f7WiU j + IB ! + 1 



Q J ui + Qi^' - 2 (t/ + y) 



-- (1 - \2U) p ti - - (aU,i + \gv 2 V, 



(52) 



(53) 



where 



a = g 



l + ^{v 2 + 2V + U+ V 



(54) 



For Qi = = Hj, V = U, a = 1 and 7^ = % these equations reduce to Eqs. (64), (67) in Chandrasekhar [5]. 
Equation (52) can be written in another form as 



A (°vy + -(<?•«%+! 



i(Q' |i + n;^) + 4| + Iv.v)n + ( 3 5 + iv.v)p 



a 1. 



+ gT- x O-^ = 0, 



(55) 



where 



* v 9 ~o 
e = £»^-^m = Q 



(56) 



with a 3 y / 7 the background part of \/—g. This corresponds to Eq. (117) in Chandrasekhar [5]; see also Eq. (44) in [6]. 
The mass conservation (continuity) equation = (gu c ). c = g + Og gives 
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= c(gu c ). c = c— = (g^f^u ) 

V — 9 > c 

= l(a 3 Q*)+-(g*v i ) +qT- 1 0~ 4 
Thus, if we assume the mass conservation, Eq. (55) gives 

The specific entropy is introduced as TdS = d (fl/c 2 ^j + p/c 2 d (l/g), thus along the flow we have 



(57) 



(58) 



TS = 



n+p (!/#■ 



_9 1 

dt a 



v v n 



+ 
a a 



1/0-2 



(59) 



This shows that for an adiabatic fluid flow the LHS of Eq. (58) vanishes; this also makes the RHS to vanish, which is 
naturally required for an ideal fluid. According to Chandrasekhar [6]: "the conservation of mass and the conservation of 
entropy are not independent requirements in the framework of general relativity. And the reason for their independence 
in the Newtonian limit is that in this limit ("c 2 — > oo") [our Eq. (57)] reduces simply to the equation of continuity." 
In [22] Blanchet, et al. suggested to use the following new combination of variable instead of Vi 



1 - c a 4 ^7 Q* 1 ' 



which becomes 



-v 2 + U + 2V + U+-) - Pi + - (Qi + IlyV) 



1 

^ =v l + — . , , 

Notice that v* is directly related to the ADM flux vector in Eq. (51) 

Ji = agev*. 

Equation TP. b — can be written as follows 



+ LT- 1 Q- A . 



(60) 



(61) 



(62) 



-9fc 



V7 



V7 



- 1 V^^ fiafc- 

,, 2 ^7 ^ 



U.i(r- I II - V... i v z + *'-) -vJPjj + y-' 



+ L 2 T -2 -4 



and we have 



-9nr 



V7 



f? = a 3 v 3 g*v* + 



V7 



p6{ + Ilj + 1 (q^ - 2VU? - lW« fc ) + L7 



(63) 



(64) 



Thus, Eq. (63) can be arranged in the following form 



1 



1 



- (av*)' + -v* u v 3 = - I U a + — 



1 



1 



Ui { -v z - U + n + V t u + 3- - v' J P, 



.2 , ,P 



1 

g*a 

1 «i 
-| - 

c 2 g* 



i + \ (sv - u)] \ P 5{ + n{ + 1 (q^ - 2vnj - rW^)] } 

C C V ' \ 

a \ '\j \ot a J \ a a J 



+ LT- z O 



2$ 



2/o-4 



(65) 



where the last line vanishes in an adiabatic ideal fluid. 
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To the O order, Eqs. (52), (53) give the Newtonian mass and momentum conservation equations, respectively 

-^(a 3 g)' + -V l (gv l ) =0, (66) 

- {avij + -v j VjVi + — h lP + Vj-n^ - - V l U = 0. (67) 
a a ag V / a 

Thus, we naturally have the Newtonian equations to OPN order, e.g., see [23]. This can be compared with the situation 
in perturbation approach where the Newtonian correspondence can be achieved only after suitable choice of different 
gauges for different variables, thus being a non-trivial result even to the linear order in perturbation approach [3,24]. 
Equations (66), (67), as well as all the 1PN equations in this section, are fully nonlinear. 

In the Friedmann background we set the PN variables U, Pi, V, Vi, Qi, and lTy equal to zero, and set g = Qb, 
II = I!;,, and p — pb, with fi = g (c 2 + n) and Ub = Qb [l + ^ (n fc + Pb/Qb)] ■ Equations (52), (57) give 



fib + 3- (fib +Pb) 



0, 



a 

gb + 3-gb 

a 



0. 



By subtracting the background part, Eq. (52) becomes 



[a 3 (<7-<J 6 )]' + 



1 



av l + 



i 

+ -^g 



V + -v l (V - U) 
a 



P-Pb 



(68) 
(69) 



gT O' 



0. 



(70) 



B. Einstein's equations 



We take Einstein's equations in the form 



rta 



8irG 



c 4 \'° 2' 
Dimensions are as follows 

\g a b] = 1, [Rob] = = [R ab ] = L- 2 , [f ab ] = [f] = [gc 2 ] = ML~ X 1 2 , 
where M indicates the dimension of mass. To 1PN order, Eqs. (8), (20) give 



A+ 4j I 3- + -zU + 4irGg 
a cr 



U- 1 (U - V) ti + 2VAU - 2A$ - (aP^^j + 8nGg fv 2 



2[V+-U 

a 



2K , , 
— - A U + - 



^ (P\ 0l - AP - 2KP?j - 8irGgavi 



0. 



a a 2 A + 4AT „ 

- + 2— — V-AwGg 

a a A a z 



[A]=L- 2 , [Gg]=T- 



-n 



2g 



5) + ^{U-V) 



-,8-KGg 



V l V4 + 



1 



n 



5 ) 



l. 



-n! 



o, 



(71) 



(72) 



(73) 
(74) 

(75) 



where Eqs. (73), (74), and (75) are R%, R°, and R) parts of Eq. (71), respectively. 

We kept the 0~ 4 term on the RHS of Eq. (75) in order to get the correct Friedmann background equation. To be 
consistent, we also need to keep the LHS to 0~ 4 which will explicitly involve 2PN order variables. But we do not 
need such efforts for our background subtraction process because all the 0~ 4 order terms involve PN order variables, 
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which do not affect the Friedmann background. To get the Friedmann background we set U, $, Pi, V, Vi, Qi, and 
riy equal to zero, and set g — gb, II = Tib, and p = pb- Then, Eqs. (73), (75) give 



Thus, we have 



a 4irG 
a 3 



ft l + T 



c z a a z 



'iPb~ 


Ac 2 ' 


c 2 _ 




3 


nA 


Pb' 


+ 


c 2 J 


c 2 _ 



a 
a 



4ttG 



3 

8ttG 



ft 1 + 



ib 



-ft i 



n, 



= o, 

2ifc 2 
^ 2_ 



Ac 2 \ 5) = 0. 



+ 



ipb 



+ 



Ac 2 



3 ' 



Kc 2 



which are the Friedmann equations. 

By subtracting the background equations, Eqs. (73)- (75) give 



~3~' 



(76) 
(77) 

(78) 
(79) 



^U + 4irG(g- Q b ) 



4 < 3V- + 3- f t> + 2V) + 6-U - 4 (U - V) . + 2VAU - 2A$ 
c 4 a V / a a 2 L '* 



aP 1 , 



+8ttG 



(gU- g b U b ) + - (p-pb) 



= 0, 



2 {^ + l u ) + ^ ( pJ i^ ~ Ap * ~ 2Xp ) ~ 8?rG ^ 



A + 4K 



V + AttG (g - g b ) 



To Q- 2 order, Eq. (80) gives 



a 2 



-4ttG (g - g 6 ) . 



(80) 

(81) 
(82) 

(83) 



which is Poisson's equation in Newton's gravity. Notice that in the Newtonian limit of our PN approximation the 
homogeneous part of density distribution is subtracted. This revised form of Poisson's equation is consistent with 
Newton's gravity, and in fact, is an improved form avoiding Jeans' swindle [25]. The decomposition of Eq. (82) into 
trace and tracefree parts gives 



V = -4ttG (q - g b ) . 

a 2 

{U-Vf i =KV8). 



Thus, compared with Eq. (83), for K = we have 



V = U, 



(84) 
(85) 

(86) 



where we ignored any surface term S with AS = and S^j = 0. Notice that we have V — U even in the presence 
of anisotropic stress; this differs from the situation in the perturbation approach where U and V in the zero-shear 
gauge (setting P 1 ^ = 0, see later) are different in the presence of the anisotropic stress [3]. For later use we take a 
divergence of Eq. (81) which gives 



A 



I* 



(87) 



Here, we reached a point to address the issue related to the background curvature in our PN approach. By taking 
a divergence of Eq. (81) and using Eqs. (83), (84), and (66) we have 
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which is apparently inconsistent. This is because the K term in Eq. (79) is related to the c -2 higher order terms in 
the background Friedmann equation. We can check this point by expanding the equation to 2PN order; in that order 
we can find that the above term can be removed by subtracting the background equation. Thus, it looks like our PN 
approximation is properly applicable only for K = 0. In this sense, our K terms do not have significance because 
these are related to (D~ 2 order higher terms through the background equation. At the moment, it is unclear whether 
this limitation of our 1PN approach to a nearly flat background is due to our subtraction process of the background 
equations (which spread over different PN orders), or intrinsic to the whole PN approximation. Meanwhile, recent 
observations favour the flat background world model with non- vanishing cosmological constant. We include the 
cosmological constant in our PN approximation which appears only in the background Friedmann equations. Ignoring 
K term we simply have V = U. 



C. ADM approach 



The ADM equations are [21,3,17] 



2 16ttG . 

—j 
3 

2 - 8ttG 

-J T,- wi v- 1 i »r:i i\r— 1 irij ir ^ 1^2 "* " 1 ' 



R W =R ij R ±k 2 + ^-E + 2K, (89) 
3 c 4 

H.i-\Ki = ^u (90) 



K^N' 1 - K ti N % N + N'-^N' 1 - K %3 K^ - -K 2 — (E + S) + A = 0, (91) 

3 c 

K^qN- 1 - K l j:k N k N- 1 + KjkN^N- 1 - K i k N k . j N~ 1 

= KK\- (Vv - + fiWJ - ^S), (92) 

EflN- 1 - EjITN- 1 - K (e + - S ij K i:j + iV" 2 (N 2 X) . . = 0, (93) 
Ji.oAT" 1 - Ji-.jWN- 1 - JjW^N- 1 - KJi + EN^N- 1 + S{ :j + SjNjN -1 = 0. (94) 

Using the ADM quantities presented in Sec. HE we can show that Eqs. (89)-(94) give the same equations we already 
have derived from Einstein's equations and the energy and momentum conservations. Equation (89) gives Eqs. (79), 
(84). Equation (90) gives Eq. (81). Equation (91) gives Eq. (73). Equation (92) gives the tracefree part of Eq. 
(82). Equation (93), (94) give Eqs. (52), (53), respectively. The ADM equations are often used in the PN approach 
[22,28,13,26]. We present the ADM approach because the ADM equations show the fully nonlinear structure of 
Einstein's gravity in a form suitable for numerical treatment. In fact, Eq. (65) can be derived from Eq. (94) using the 
identification made in Eq. (62). 



IV. 1PN EQUATIONS 
A. Complete equations to 1PN order 



Here we summarize the complete set of equations valid to 1PN order in the cosmological situation. The background 
variables a, Qf,, IIj and pb are provided by solving Eqs. (68), (69), (78), and (79) 



a 


4ttG 


a 


3 


a 2 

a 2 ~ 


8?rG 
3 Qb 



Qb 1 



lb, 

A 



3pb 
n 1 



Ac 2 
~3~' 



IT, 



fi b + 3- (pb+Pb) = 0, 

Qb + i-Qb = 0, 

a 



Kc 2 Ac 2 

— o- + 



3 ' 



(95) 
(96) 

(97) 
(98) 
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where = QbC 2 (l + Ilb/c 2 ). In Eqs. (95)-(97) only two are independent, and from Eqs. (97), (98) we have Qbtlb + 
3(d/a)p b = 0. 

To the Newtonian order we have Eqs. (66), (67), and (83) 



-(a 3 e ) + -v< (ev i ) = o, 



- {avi)' + -v j VjVi + — (v iP + VjU{) - -ViU = 0, 
a a ag \ ) a 



a 
A 



U + AttG (g - g b ) = 0. 



(99) 
(100) 
(101) 



In the Newtonian case energy conservation and mass conservation provide an additional equation which is Eq. (58) 

(I + r ■ v ) n + ( 3 ; 4 V ■ ") ; + 5 («• h + n H<) = »■ <™> 

In the non-expanding background this equation gives the well known energy conservation equation in the Newtonian 
theory, see for example Eq. (2.36) in [27]. The Newtonian order hydrodynamic equations are valid in the presence of 
general background curvature K . 

In the following we set K — to the 1PN order, because our PN expansion fails to apply in the presence of K term. 
Thus we have V — U . The energy and the momentum conservation equations valid to 1PN order are derived in Eqs. 
(70), (53), and these are 



(aVT + 



1 



(103) 



cfVi + — (Qi + UtjV 1 



+ - 



a 
1 

c z 



;2U p, 



1 

+ - 

a 

1 



iV j + U{ + 1 (cp Vi + QiV j - AUufj 



a + —gv 2 U, 



^{§- t + l--v)u + l (u 2 - *) . - 1 {aPi y - 2 -^p m + ±u d ui 



where 



a = g 



v 2 + 2U + U+- 



(104) 



(105) 



To the Newtonian order Eqs. (103), (104) for g and Vi reduce to Eqs. (99), (100), respectively. From Eq. (55), (65) 
we can derive alternative forms 



- (a 3 ,*)' + - (q*v% = -- 2 [- (V |j + ny „) + g 

U.i + - 



d_ i 

dt 



v-V)II+(3- + -V-v)p 



(106) 



- (av*)' + -v* u v J = - 



\v 2 - U + n + 3^ ) U,i + 2*,* - yiPju 



+ 



a \ k J \j \dt a J \ a a 



c 2 g* 



(107) 



where 



g = g 



1 + 



1 (I 



-v 2 +3£7 



■U,; = Vi + 



\ v 2 + 3U + n + Vi - p l + - (Qi + iii 



(108) 
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The hydrodynamic and thermodynamic variables II, p, Qi and Ily should be provided by specifying the equations of 
state and the thermodynamic state of the system under consideration. 

The metric variables U, <E> and Pi can be expressed in terms of the Newtonian variables g, Vi, IT, p, Qi and 11^ by 
using Einstein's equations. Equations (83), (84) give 



-^U = -4nG(g-g b ), 



(109) 



which already determines U. Equations (80), (81), and (87) give 

2A$-2£/AC/ + (aP\ t 



A 11 
-U + 4*G(Q-e b ) + -z\- 5 



3U + 9-U + 6-U 
a a 



+8irG 



gv 2 + -( e U- g b n b ) + ^{P- Pb) 



0. 



A 

a 



= 4ttG 

a z 
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a \a 
a 



—P l = -WirGgvi + - ( -P J |? +AU + 4-U 



- {qv 1 ) h + -{Q- Qb) 



(110) 
(111) 
(112) 



Notice that to 1PN order Einstein's equations do not involve the anisotropic stress or flux term. To 1PN order Eq. 
(110) determines $. Our conservation equations (103), (104) contain U terms. In order to handle these terms it 
was suggested in [28] that the Poisson-type equation in Eq. (112) provides better numerical accuracy. Notice that, 
whereas only the spatial gradient of the potential U appears in the Newtonian limit, we have bare U terms present 
in the 1PN order [29]. 

In order to handle these equations we have the freedom to take one temporal gauge condition. This corresponds 
to imposing a condition on P 1 ^. or $. Gauge related issues will be addressed in detail in Sec. V. There, we will show 
that all our variables in this section are spatially gauge-invariant, and are temporally gauge-ready. 



B. Ideal fluid case 



We consider an ideal fluid, i.e., Qi = = 11^ , and assume the adiabatic condition Eq. (102) applies. Thus, the 
internal energy is determined by the energy conservation equation 



at a J \ a a J g 



(113) 



To the background order we have Eqs. (95)-(98). To the Newtonian order we have Eqs. (99)-(102). Equations (103), 
(104) become 



1 ' 3_v , 1 , 1 „ (tt , a ,.2 P 



(a^y + ( av >) + g t U+ v 



\ (a^v,)' + - {avy), - - [c7 + \gv 2 ) U. l + - ( 1 - X{2JJ ) p,, 
a 4 v ' a y ' U a \ c z J a \ c z ) 



1 



d 1 



= 0, 



where a is given in Eq. (105). These equations can be written as 



- (avi) + -Viuv J U i H h 

a a u a a g c l 



X - V 2 + 4U + U + P)^ 

a V QJ Q 



(114) 



(115) 



(116) 



= 0. 



(117) 
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Alternative forms follow from Eqs. (106), (107) 



-l(aV) + ^V) 



0. 



a 



1 + 



1 /3 



where 



= Q 



1 + 



1 /l 



-v 2 - u + u + 



i 



1 1 



(2*,i-^|i) 



-u 2 + 3J7 + n + i- )vi - Pi 



(118) 
(119) 

(120) 



All these equivalent three sets of equations are written in Eulerian forms. Einstein's equations in (110)-(112) will 
provide the metric perturbation variables U, $ and Pi in terms of the Newtonian fluid variables (g, v l , II, and p), 
after taking one temporal gauge condition. The pressure term should be provided by an equation of state. In the 
case of a zero-pressure fluid with vanishing internal energy we can set p = = II. We can take any set of equations 
of motion depending on the mathematical convenience in numerical treatments. 



V. GAUGE ISSUE 



A. Gauge transformation 

We consider the following transformation between two coordinates x a and x a 

x a = x a +C(x e ). 

For a tensor quantity we use the tensor transformation property between x a and x a coordinates 



tab(x e ) 



dx" dx b cd(x h 



(121) 



(122) 



Comparing tensor quantities at the same spacetime point, x a , we can derive the gauge transformation property of a 
tensor quantity, see Eq. (226) in [17], 



t ab (x e ) = t ab (x e ) - 2t c{a e M - t ab j c + 2t c(a i d M e. d + t c d c J\ b 

( 2i C ,( Jb)c,d + 2tc(ai C ,b)d + 7;tab,cd£ C + tab,c£ C , d 



(123) 



We considered £ a to the second perturbational order which will turn out to be sufficient for the 1PN order. In the 
following we will consider the gauge transformation properties of the metric and the energy-momentum variables. 
As the metric we consider the following more generalized form 



.9oo = - 
9oi = -■ 



9ij = a 2 



1 



1 



-2U- 



1 



aPi + 0~ 

1 - i^" ; 



(2C. 



o- 



2C W) + 2C K 



(124) 



where Cj is transverse (C 4 ^ = 0), and is transverse and tracefree {Cl - = = Q); indices of C, and Cjj are 
based on 7^. In Eq. (124) we introduced 10 independent metric components: U and $ together (1-componcnt), Pj 
(3-components), V (1-component), C (1-component), Cj (2-components), and CV, (2-components). It is known that 
gravitational waves show up in the 2.5PN order [8]. Thus, we ignore the transverse and tracefree part, i.e., set CV, = 
to 1PN order. 

We wish to keep the metric in the form of Eq. (124) in any coordinate system. Thus, we take the transformation 
variable £ a to be a PN-order quantity. We consider coordinate transformations which satisfy 
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c e 3 cr a 

where the indices of the £( 2 )"s are based on 7^ . 

The gauge transformation of Eq. (123) applied to the metric gives 

u = u + i^°, 

$ = $ _|_ 1^(4)0 _ 1^(2)0 _ .^(2)i _ 1^(2)0^(2)0 _ 1^(2)0^(2)0 _ 1 / l^(2)i^(2)0 

2 2 2ct 2 4 2 V ci 



(125) 



,(2)0 



^ = J p i --e (4)0 i + a(-e 



t(2) 



jj + M2)o\ t(2)0 + 1^(2)0 ,(2)0 + J_ / t(2)jt(2)0 \ 



(A + SK) AC = (A + 3J0 (AC - + ^ (^ (a)0 , 1 ) - ^ (W) 



i A f 2/v )C = ( A 4- 2 A') | a - ) - — ( 



J_t(2)i 
3a l> r, 

1 



+ i Vi A-{|(A + 3iO€% + i 
(A + 3A") y = (A + 3K) (v - -C (2)0 



(2)j 

a 

A + 2K 

4a 2 



.(2)0 



^0,^)0 .) + _1_ ( ? (2)0 /( 2)0^ 



(126) 
(127) 
(128) 
(129) 

(130) 



(131) 
(132) 



Equations (126), (127) follow from the transformation of g 00 ; Eqs. (128), (129) follow from g 0l ; Eqs. (130)-(132) follow 
from gij . 

Equation (128) shows that £( 2 ) is spatially constant, 



;(2)o _ t(2)0 



After this, Eq. (130) gives 



£™ u (t). 



AC = AC - -£ (2)i ,. 



(133) 



(134) 



Thus, by choosing C = (AC = is enough) as the gauge condition (this implies that we set C — to be valid in 
any coordinate), we have 



(2)i _ 



= 0. 



(135) 



Imposing the conditions in Eqs. (133), (135), Eq. (131) gives 

1 



(136) 

Thus, by choosing C, = to be valid in any coordinate, i.e., choosing d = as the gauge condition, we have 



d 2) = o. 

Under the conditions in Eqs. (133), (137) the remaining metric perturbation variables transform as 

U = U + £ (2)0 , V = V - -£ (2)0 , $ = $ + -C (4)0 - -U£ (2)0 - i^ 2 )°^2)° - 1^(2)0^(2)0, 

ci 2 2 2 4 

p _ p _ 1 t(4)0 



(137) 



(138) 



Notice that for non- vanishing £( 2 )°, U and U transform differently under the gauge transformation even in the flat 
cosmological background. Since the spatially constant £( 2 )° can be absorbed by a global redefinition of the time 
coordinate, without losing generality we can set £( 2 )° equal to zero 
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£ (2)0 = o, 



(139) 



thus allowing V = U in any coordinate. If we set ^ 2 ^° = but do not take the spatial gauge which fixes £- 2 ^ , from 
Eqs. (126)-(132), the metric variables transform as 

U = U, (140) 

3> = $ + k (4) °-f U,^\ (141) 
I La 

A = fl-^ W °,i + «(^ (a) )', (142) 

AC = AC - -J 2) \ v (143) 

(A + 2K) Ci = (A + 2K) (Ci ^f) - + ^A" 1 (A + 3X)C ( %, (144) 

V = V. (145) 

Taking the spatial gauge conditions C = = Ci we have £,- 2 ^ = = an d the remaining metric perturbation 

variables transform as 

U = U, V = V, * = * + U W °, Pi = Pi- -& )0 ,i. (146) 
In the perturbation analysis we call C = = Ci the spatial C-gauge [17]. Apparently, under these gauge conditions 

(2) (2) 

the spatial gauge transformation function to 1PN order, Q , is fixed completely, i.e., Q =0. By taking such 
gauge conditions, the only remaining gauge transformation function to 1PN order is £( 4 ) . This temporal gauge 
transformation function affects only <!> and Pi. If we take $ = as the temporal gauge condition, we have £( 4 )° = 0, 
thus £( 4 )° does not vanish even after imposing the temporal gauge condition and has general dependence on spatial 
coordinate £( 4 )°(x). Thus, in this gauge, even after imposing the temporal gauge condition the temporal gauge mode 
is not fixed completely. Whereas, if we take P 1 ^ = as the temporal gauge condition, we have £ (4)0 = 0, thus the 
temporal gauge condition is fixed completely. In fact, from the gauge transformation properties of $ and Pi we can 
make the following combinations 

*,i + l(aPiy, A$ + i(aP* K )', (147) 

which are invariant under the temporal gauge transformation, i.e., temporally gauge- invariant. Meanwhile, U and V 
are already temporally gauge- invariant. Since our spatial C-gauge also has removed the spatial gauge transformation 
function completely, we can correspond each remaining variable to a unique (spatially and temporally) gauge-invariant 
combination. Thus, in this sense, we can equivalently regard all our remaining variables as gauge-invariant ones. For 
example, for K = 0, from Eqs. (140)-(145) we can show 



Pi 



a(di + C,ij ' =Pi + a{Ci + C,i)' - ^ (4) °,»- (148) 



Thus, Pi + a (Ci + C.i)' is spatially gauge-invariant and becomes Pi under the spatial C-gauge. This implies that 
Pi under the spatial C-gauge is equivalent to a unique gauge- invariant combination Pi + a (C + Cj)'. Similarly, 
in the case of the temporal gauge, from Eq. (147) we can show that A$ under the P 1 ^ = gauge is the same as 

a unique gauge-invariant combination A<I> + |(aPV)'. In this sense, under gauge conditions which fix the gauge 
mode completely, the remaining variables can be regarded as equivalently gauge-invariant ones with corresponding 
gauge-invariant variables. 

In Eq. (1) we began our 1PN analysis by choosing the spatial C-gauge 

C = = C (149) 

By examining Eqs. (140)-(145) we notice that the spatial C-gauge is most economic in fixing the spatial gauge mode 
completely without due alternative. In this sense the spatial C-gauge can be regarded as a unique choice and we do 
not lose any mathematical convenience by taking this spatial gauge condition. In the literature these conditions are 
often expressed in the following forms. The spatial component of the harmonic gauge condition sets the 1PN part of 
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-abfi 



ab 



99 lc 



= \^{C- 1 + C t )+L- 1 0-\ 



(150) 



equal to zero. We can also set the 1PN part of 



~9° k (§ij,k - \~9jk^ = (ic ti + C^j + L- l O-\ 



(151) 



equal to zero; notice that in Eqs. (150), (151) we assumed K = 0. In either case we arrive at Eq. (149). We have 
shown that under the conditions in Eq. (149), the spatial gauge transformation is fixed completely without losing any 
generality or convenience. Whereas, we have not chosen the temporal gauge condition which can be best achieved by 
imposing a condition on P 1 ^. We call this the gauge-ready strategy [18,17]. It is convenient to choose this remaining 
temporal gauge condition depending on the problem we encounter; or to try many different ones in order to find out 
the best suitable one or ones. 

Now, let us consider the gauge transformation property of the energy-momentum tensor. We set ^ 2 ^° = 0. The 
gauge transformation property of Eq. (123) applied to the energy-momentum tensor in Eq. (20) gives 



1 



(2) 



a q g 
1 



Vi = 



Qi = 



(2) 



v[ 2) - a i -e 



(2) 



c z a \ 3 



(2)i 



-nv (2)j 



1 1 

c 2 a 



^(^)-^^ (2 V)+ n .l^ (2)fe + 2n M^ (2 %) 



(152) 
(153) 
(154) 
(155) 
(156) 
(157) 



Equations (152), (153) follow from the transformation property of T 00 ; Eqs. (154), (155) follow from T oi : Eqs. (156), 
(157) follow from Ty. The transformation functions and are not determined, see below. 
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If we take the spatial C-gauge, i.e., set Q ' = 0, we have 



q = q+\q {2 \ 



n = n 



,(2) 



1 (2) 



(2) 



P = p + o~ 



(158) 



Thus, we have 



1 + ^ 



q U + ^n 



ii- 11 
-o-Qi = Vi + -^-Qf 

&■ g c £ g 



(159) 



Although the gauge transformation properties of g and II are not determined individually, the energy density 
[i = £>(l + n/c 2 ) is gauge-invariant under our C-gauge. For vanishing II we have g^ = 0. Similarly, the gauge 
transformation properties of Vi and Qi are not determined individually. For vanishing flux term in any coordinate, 
Qi = 0, we have = 0, thus v t = Vi + <D~ A . Thus, the gauge transformation property of Vi is determined only for 
vanishing flux term. From Eq. (108) we have 



„. + i(W^„, 

c z \a ' g 



(160) 



We can check that up to the 1PN order Einstein's equations and the energy and momentum conservation equations 
in Sec. IV are invariant under the gauge transformation. In the following we introduce several temporal gauge 
conditions each of which fixes the temporal gauge mode completely. In relativistic gravity, gauge conditions consist of 
four non-tensorial relations imposed on the metric tensor or the energy-momentum tensor. Our purpose is to employ 
the temporal gauge (slicing) condition to make the resulting equations simple for mathematical/numerical treatment. 
Depending on the situation we can also take alternative spatial gauge conditions for that purpose. In the following 
we impose the spatial C-gauge in Eq. (149), and set K = 0. 
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B. Chandrasekhar's gauge 



Chandrasekhar's temporal gauge condition in his Eq. (24) of [5] corresponds to setting the 1PN part of 



9 V [goi,j ~ \ 



c a 



1 (1 



equal to zero; in the literature this is often called the standard PN gauge [22] . We take 

i 

\i 



-P i H + 3U + m-U = 0, 
a 1 a 



(161) 



(162) 



as Chandrasekhar's gauge. We used the freedom to add an arbitrary m{a/a)U term with m, a real number. In this 
case Eqs. (110), (111) give 



4« 



-WnGgvi + — 
a 



U-(m-A)-U 



(163) 



A 1 

— U + 4nG(g- g b ) + — 



2^*-(m-3)-U + 



(a \ a a 
(6 — to) TO — 



U 



+8nG 



qv 2 + x - (gH - g b n b ) + U(g- Q b ) + \{p- Pb ) 



0. 



(164) 



Therefore, U, Pi, and <!> are determined by Eqs. (109), (163), and (164). The variable U can be determined from Eq. 
(112). This completes our 1PN scheme based on Chandrasekhar's gauge. When we handle this complete set of 1PN 
equations numerically, we should monitor whether the chosen gauge condition is satisfied always; this could be used 
to control the numerical accuracy. 



C. Uniform-expansion gauge 

The expansion scalar of the normal-frame vector, 6 = n c . c , is given in Eq. (36). It is the same as the trace of 
extrinsic curvature K with a minus sign, see Eq. (50). Taking the 1PN part of K equal to zero 

-PV + 3J7 + 3-f/ = 0, (165) 
a 1 a 

can be naturally called the uniform-expansion gauge. In the literature it is often called the ADM gauge [22], or the 
maximal slicing condition in numerical relativity [30]. This condition corresponds to the to = 3 case of Chandrasekhar's 
gauge in Eq. (162). 



D. Transverse-shear gauge 



The shear of the normal-frame vector is given in Eq. (37). Thus, the gauge condition 



o. 



can be called the transverse-shear gauge. In this case Eqs. (110), (111) give 



4« 



-WnGgvi + -{U+-U 



A 



-fU + 4ttG (q - Q b ) + ^<2—^> + 3U + 9-U + 6-U 



(166) 
(167) 



+8ttG 



gv 2 + -(qIL- g b Tlb) + U (g - g b ) + - (p - p b ) 



0. 



(168) 



Therefore, U, Pi, and <f> are determined by Eqs. (109), (167), and (168). U can be determined from Eq. (112). We 
still have U which can be determined similarly by taking the time derivative of Eq. (112) and by using the background 
and the Newtonian equations in Eqs. (95)-(101). 
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E. Harmonic gauge 



We have 



1 d 
-3- 

c a 



1 

^3 



1 



-F' 



-4(7- 



■6 a -U 
a 



+ L- 1 0- 5 . 



-g v ' / , a c a c° \ a 

The well known harmonic gauge condition sets the 1PN part of Eq. (169) equal to zero, thus we take 

1 



AU- 



m—U : 
a 



0. 



(169) 



(170) 



where we used a freedom to add an arbitrary m(a/a)U term. Combined with Eq. (150) the full harmonic gauge 
condition can be expressed by setting the 1PN parts of g ab T c ab equal to zero. In the perturbation approach the 
harmonic gauge (often called de Donder gauge) is known to be a bad choice because the condition involves derivatives 
of the metric variables; this leads to an incomplete gauge fixing and consequently we would have to handle higher 
order differential equations which is unnecessary, see the Appendix in [31]. In this gauge Eqs. (110), (111) give 



A 7 P t = -16irGgv t - (m - 4) --U ti , 



*U + l*G{e-Q b )+U2±*. 



a a 



(171) 



U-(m- 



i) a -u + 

a 



(6 — m) m- 



U 



+8nG 



gv 2 + - (qIL - g b n b ) + U (g - g b ) + - (p - p b ) 



0. 



(172) 



Therefore, U, Pi, and <& are determined by Eqs. (109), (171), and (172). The variable U can be determined from Eq. 
(112). We also have a U term in Eq. (172) which might demand for a more involved numerical implementation as 
explained below Eq. (168). However, its presence makes the propagating nature of the 1PN order metric fluctuations 
apparent in this gauge condition. This time-delayed propagation of the gravitational field in the 1PN approximation 
can be compared with the action-at-a-distance nature of Poisson's equation in the Newtonian order in Eq. (172). 
We anticipate that the relativistic time-delayed propagation could lead to a secular (time cumulative) effect. This 
could be important even in the case in which the 1PN correction terms are small compared with the Newtonian 
terms; in the large-scale clustering regions we would expect (v/c) 2 ~ GM/(Rc 2 ) ~ 10~ 6 ~ 10~ 4 , see Sec. VI. Using 
goo = — 1 + 2U/c 2 we have 



1 



2$ 



i - 



(173) 



thus 



□U = U ;C C = — U 
a 2 

= 4t/ 



and Eq. (172) can be written as 



1 



□U + AttG (g - g b ) 
1 fA 



u + 3-u + 2u4u^l + r- 2 e>- 4 

a a 2 J 

4 (U 2 - 2$) + U + 3-U + 2U^U 



+ T~ 2 Q- 



(174) 



+ 3 72 U "("»■ 



4)^U- 
a 



a a 

(6 — m) m^r 

a a z 



U + 8ttG 



gv 2 + ^(qU- g b Il b ) + ^ (p - p b ) 



0. (175) 



Thus, in the harmonic gauge condition the propagation speed of the gravitational potential U is the same as the 
speed of light c. (When we mention the propagation speed, we are considering the D'Alembertian part of the wave 
equation ignoring the nonlinear terms and background expansion. In this case we have DU = a~ 2 AU — c~ 2 U, and 
for K = 0, A becomes the ordinary Laplacian in flat space.) Apparently, the wave speed of the metric (potential) can 
take an arbitrary value depending on the temporal gauge condition we choose. For example, if we take 



-PY- +nU + m-U = 0, 



(176) 
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as the gauge condition, with n and m real numbers, Eqs. (110), (111) give 



4« 

a 2 



-WnGgvi 

a 



(n-4)t/ + {m-A)-U 



^rU + 4ttG - g b ) + \ \ 2^$ - (n - 3) £/ - (2n + m - 9) -U + 



(6 — TO) TO— - 



+87rG 



e« 2 + - ( e n - gb n b ) + u( g - g b ) + -(p- Pb ) 



0. 



(177) 



(178) 



In this case, for n > 3, the speed of propagation corresponds to 

c 



\/n — 3' 



(179) 



which can take an arbitrary value depending on our choice of the value of n. It becomes c for n = 4 (e.g., the harmonic 
gauge), and infinity for n — 3 (e.g., Chandrasekhar's gauge and the uniform-expansion gauge). In the case of the 
transverse-shear gauge we have n = 0, thus Eq. (168) is no longer a wave equation. 



F. Transformation between two gauges 

Now, let us show how we can relate the equations and solutions known in one gauge condition to the ones in any 
other gauge condition. Our spatial C-gauge condition already fixed the gauge transformation function Q ' = 0, and 
we have £ (2)0 ee 0. Under the remaining (temporal) gauge transformation 

t = i+lc (4) °, x*=x\ (180) 

we have 

$ = $ + Ie (4)0 , P 4 = P,--£ (4)0 4 , (181) 
A a 

and all the other variables are gauge- invariant. As an example, let us consider the two gauge conditions used in Sec. 
VD and Sec. VE. Let us assume that the x a coordinate is the transverse-shear gauge in Sec. VD, and x a is the 
harmonic gauge in Sec. VE. Thus, we have P 1 ^ ee in x a , and 4{7 + to|£7 + a~ 1 P t ^ i ee in x a . From Eq. (181) we 
have 

W + m-u] ={)- 4e (4) °, (182) 



a J a 2 



thus 



4C (4) ° = 4f7 + m-U. (183) 



a 2 ' a 



Thus, when we transform the result known in the transverse-shear gauge to the harmonic gauge, we use Eq. (181) 
with £( 4 ) given in Eq. (183). All the other variables are invariant under the gauge transformation. We can show that 
under this transformation Eqs. (167), (168) become Eqs. (171), (172), respectively. In this way, if we know the solution 
in any gauge condition the rest of the solutions in other gauge conditions can be simply derived without solving the 
equations again. Thus, as in other gauge theories (like Maxwell's or Yang-Mills' theories) the gauge condition should 
be deployed to our advantage in handling the problem mathematically. 



VI. DISCUSSION 



Here we compare the PN approach with the pcrturbative one. The perturbation analysis is based on the perturbation 
expansion of the metric and energy-momentum variables in a given background. All perturbation variables are assumed 
to be small. In linear order perturbations we keep only the first-order deviations from the background [1-3], whereas 
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in the weakly nonlinear perturbations we keep higher-order deviations to the desired order [17]. Contrary to the PN 
approach the perturbation analysis is applicable in the strong gravity regime and on all cosmological scales as long 
as the perturbations are linear or weakly nonlinear. 

Meanwhile, in the PN approach, by assuming weak gravitational fields and slow motions, we try to provide the 
general relativistic correction terms for the Newtonian equations of motion. Thus, in the PN approach, in fact, 
we abandon the geometric spirit of general relativity and recover the concept of absolute space and absolute time. 
Although this could be regarded as a shortcoming of the PN approach, in this way it provides the relativistic effects 
in forms of the correction terms to the well known Newtonian equations, thus enabling us to use simpler conventional 
(numerical) treatment. We expand the metric and the energy- momentum variables in powers of v/c in a given 
background spacetime. In nearly virialized system we have GM/(Rc 2 ) ~ (v/c) 2 which is assumed to be small. 
Thus, no strong gravity situation is allowed and the results are valid inside horizon ^GM/(Rc 2 ) ~ R/(c/H) < 1. In 
comparison to the perturbation approach, however, the resulting equations in the 1PN approximation can be regarded 
as fully nonlinear. As in the perturbative case, even in our cosmological PN approach we assume the presence of a 
Robertson- Walker cosmological background. 

As we have summarized in the introduction, our studies of the weakly nonlinear regime of a zero-pressure cos- 
mological medium showed that the Newtonian equations are quite successful even near the horizon scale where the 
fluctuations are supposed to be near linear stage. We have shown that to the second order in perturbations, except 
for the presence of gravitational waves, the relativistic equations coincide exactly with the Newtonian ones. The pure 
relativistic correction terms appearing in the third order are independent of the presence of the horizon and are small 
compared with the second-order terms by a factor ST/T ~ 10 -5 , thus negligible. 

In order to properly estimate the relativistic effects in the evolution of large-scale cosmic structures we have to 
implement our equations in a hydrodynamic cosmological numerical simulation. The PN correction terms are (v/c) 2 
or GM/(Rc 2 ) orders smaller than the Newtonian terms. In Newtonian numerical simulations the maximum large 
scale velocity field of a cluster flow reaches nearly 3000km/ sec and the typical value for the velocity is about an order 
of magnitude smaller than this [32]. Thus, we may estimate the 1PN effect to be of the order (v/c) 2 ~ 10~ 6 ~ 10~ 4 , 
thus quite small. We already mentioned that the PN approximation is applicable inside the horizon only. Considering 
the action-at-a-distance nature of the Newtonian gravity theory it is important to check the domain of validity of 
the Newtonian theory in the nonlinear evolution of cosmological structures. The PN approach provides a way to find 
out the relativistic effects in such a regime. We anticipate that the propagating nature of the gravitational field with 
finite speed in relativistic gravity theory, compared with the instantaneous propagation in Newton's theory, could 
lead to accumulative (or secular) relativistic effects. In order to estimate relativistic effects it would be appropriate 
to consider a single cold dark matter component as a zero-pressure fluid without internal energy. In such a case our 
equations in Sec. IV B with p = = II and the metric variables U, $, and Pi presented in various gauge conditions 
in Sec. VB-VE would provide a complete set of equations expressed in various forms. 

Large-scale structures are still near the linear regime, and due to the enhanced amplitude of the initial mass power 
spectrum in the small scale the gravitational evolution causes nonlinear regions to begin in the small scale and to 
propagate to larger scales [33]. The current belief is that in small scale structures, where the nonlinearity is important, 
the dynamical time-scale is much longer than the light travel time over this scale, thus the time dilation effect from 
relativistic gravity is not important. The 1PN order relativistic effects could be important in the tidal interactions 
among clusters of galaxies, where the dynamical time-scale could become substantial compared with the light crossing 
time of the scales involved [29] . 

We can perform numerical simulations based on any two different gauge conditions. The results should coincide 
after making a gauge transformation between the two gauges as in Sec. VF. Also, the two simulations should give 
identical result for any given gauge-invariant variable. These might provide a way to check the numerical accuracy of 
the simulations. Analyses based on different gauge conditions would lead to the different results; this is in the sense 
that a given variable evaluated in two different gauges are actually two different variables, unless the original variable 
is gauge- invariant. Even for the temporal gauge condition (after fixing the spatial C-gauge) there are infinitely 
many different gauge conditions available. Since each of the gauge conditions displayed in Sec. VB-VE fixes the 
temporal gauge mode completely (and the spatial gauge modes are already completely fixed by our spatial C-gauge), 
all remaining variables are equivalently gauge- invariant. Thus, apparently, the gauge- invariance does not guarantee us 
to associate the variables with physically measurable quantities. The identification of physically measurable quantities 
out of infinitely many gauge-invariant candidates is still an open issue which remains to be addressed. The gauge 
invariance assures that the value should not depend on the gauge we take and this fact can be used to check and 
control the numerical accuracy of the simulations. 

In a related context, we have shown that the propagation speed of the gravitational potential depends on the 
temporal gauge condition we take, see Eq. (179). A similar situation occurs in classical electrodynamics. The 
propagation speeds of the electromagnetic scalar and vector potentials are the speed of light c in the Lorenz gauge, 
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whereas that of scalar potential becomes infinite in the Coulomb gauge. In electromagnctism the issue is well resolved 
by showing that the electric and magnetic fields propagate with c independently of the adopted gauge condition [34] . 
In our 1PN approach, however, we have not been able to resolve the case. In the PN case this is related to identifying 
the physically relevant gauge condition out of infinitely many different gauge conditions available (which are all gauge- 
invariant), an issue we have described in the previous paragraph. One difference compared to electromagnetism is that 
the PN approach addresses fully nonlinear gravitational dynamics whereas electromagnetics concerns linear processes. 
We wish to address these important issues in a future occasion. 

Just like the solar system tests of Einstein's gravity theory, the nonlinear evolution of the large-scale cosmic structure 
could provide another regime where the gravitational field is weak and the motions are slow so that the post-Newtonian 
approximation would be practically adequate to describe the ever-present relativistic effects using Newtonian-like 
equations. The problem is whether such effects are significant enough to be detected in future observations and 
numerical experiments. This important issue is left for future studies. We hope that our set of 1PN order equations 
and our strategy for using those equations would be useful for such studies in the cosmological context. In more 
realistic cosmological situations we have dust and cold dark matter which can be approximated by two zero-pressure 
ideal fluids. We can also include pressure and dissipation effects in the case of dust. The multi-component situation 
of cosmological IPN hydrodynamics will be considered in later extensions of this work. Extending our formulation to 
higher order PN approximation is also an apparent next step which would be tedious but straightforward. 
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